#ifndef GAUSS_H
#define GAUSS_H

#include <qmath.h>
#include <stdlib.h>

#define a(i,j) a[n*(i)+(j)]
#define swap(x,y,d) {d=x;x=y;y=d;}

void gauss(int,double*,double*);

/*

59.895571;30.082165
59.895571;30.253827
59.809358;30.253827
59.809358;30.082165

Xi=a11xi+a12yi+a13-a31xiXi-a32yiXi
Yi=a21xi+a22yi+a23-a31xiYi-a32yiYi

x1 y1 1 0 0 0 -x1X1 -y1X1 a11 X1
0 0 0 x1 y1 1 -x1Y1 -y1Y1 a12 Y1
x2 y2 1 0 0 0 -x2X2 -y2X2 a13 X2
0 0 0 x2 y2 1 -x2Y2 -y2Y2 a21 Y2
x3 y3 1 0 0 0 -x3X3 -y3X3 a22 X3
0 0 0 x3 y3 1 -x3Y3 -y3Y3 a23 Y3
x4 y4 1 0 0 0 -x4X4 -y4X4 a31 X4
0 0 0 x4 y4 1 -x4Y4 -y4Y4 a32 Y4

*/


/*
int main()
{
    int n=8,i;
    double a[]={
                    0,0,1,0,0,0,0,0,
                    511, 0, 1, 0, 0, 0,-308644, 0,
                    0, 255, 1, 0, 0, 0, 0,-22950,
                    511,255,1,0,0,0,-272363,-135915,
                    0,0,0,0,0,1,0,0,
                    0,0,0,511,0,1,-133371,0,
                    0,0,0,0,255,1,0,-79050,
                    0,0,0,511,255,1,-228417,-113985
                };
    double b[]={126,604,90,533,116,261,310,447};
    gauss(n,a,b);
    for(i=0;i<n;i++) printf("%.18le\n",b[i]);
    return 0;
}
*/
void gauss(int n,double *a,double *b)
{
    int i,j,k;
    double z;
    for(i=0;i<n-1;i++)
    {
        k=i;
        z=fabs(a(i,i));
        for(j=i+1;j<n;j++)
        {
            if(fabs(a(j,i))>z)
            {
                k=j;
                z=fabs(a(j,i));
            }
        }
        if(k!=i)
        {
            for(j=i;j<n;j++)
                swap(a(i,j),a(k,j),z);
            swap(b[i],b[k],z);
        }
        for(j=i+1;j<n;j++)
        {
            z=a(j,i)/a(i,i);
            for(k=i+1;k<n;k++)
                a(j,k)-=z*a(i,k);
            b[j]-=z*b[i];
        }
    }
    for(i=n-1;i>=0;i--)
    {
        for(j=i+1;j<n;j++)
            b[i]-=b[j]*a(i,j);
        b[i]/=a(i,i);
    }
}

#endif // GAUSS_H
